# -*- coding: utf-8 -*-
"""
Created on Mon Nov 22 12:56:29 2021

@author: zhiyinan
"""


import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
import Parameters_1D as par_RNN 
from Polar1D_vectorized import *
import van_Genuchten1D_vectorized as vg 
from scipy import integrate
from sklearn.preprocessing import MinMaxScaler
from casadi import *
import matplotlib.font_manager as font_manager
import time
import pandas as pd
from sklearn.linear_model import Ridge
import random
# %% Input Signal generation
np.random.seed(5)
random.seed(5)


soilPars = loamySoil() # The soil parameters
totalDepth, axialNodes, axialDistance, totalNodes, numOfActuators, interPoints = par_RNN.spatialVariables_1D()  #Spatial parameters
samplingTime, samplingTimeInternal, internalTimeSteps, totTimeSteps = \
                                 par_RNN.temporalVariable_dataGen(30000*2*60*60) #time parameters

# lbu = -5.10e-05  #  unit [m/s]
# ubu = -4.94e-05  #  unit [m/s]

# irr_freq = int(20*60*60/samplingTime)   # PRBS change once only if a minimum of 100 hr is reached

u_train = np.zeros((totTimeSteps,1))

i=0
while i < 30000:
    irr_freq = int(30*60*60/samplingTime)
    ulist = np.array([-7e-6, -6e-6, -5e-6, -4e-6, -3e-6, -2e-6, -1.5e-6, -1e-6, -0.5e-6])
    # ulist = [-2/1000/3600, -1.5/1000/3600, -0.9/1000/3600, -1.2/1000/3600, -0.5/1000/3600]
    t = irr_freq*random.randrange(1,3) + random.randrange(irr_freq)
    u_train[i] = random.choice(ulist)
    i = i+t
    print(i)

plt.figure(1)
plt.plot(u_train[:1000])                


initCondition = -0.9*np.ones(totalNodes)
#Arrays to store the results




def ODE(t, head, irrigAmount, cropCoeff, refEvap, soilPars):
    return Richards1D_vectorized(head, irrigAmount, cropCoeff, refEvap, soilPars)

cropCoeff = 0.88*np.ones(totTimeSteps)
refEvap = 3.102e-08*np.ones(totTimeSteps)

 # %% Integration 
# u_train = np.loadtxt('u_NN_train_crop_2hr_impulse_f30.txt')/50
# u_train = (np.loadtxt('u_NN_train_crop_2hr_impulse_f30_small.txt'))[:30000]
totTimeSteps = u_train.shape[0]
headArray_1= np.zeros((totTimeSteps+1, totalNodes)) # Pressure head, the state
headArray_1[0,:] = initCondition

volMoisture_1 = np.zeros((totTimeSteps+1, totalNodes)) # Volumetric moisture, the measured output
volMoisture_1[0,:] = volMoistureAllNodes_1D(initCondition, soilPars).ravel()


i=1
temp = headArray_1[i-1]
ui = u_train[i-1]
for t in range(int(samplingTime/samplingTimeInternal)): 
    sol = integrate.solve_ivp(ODE, args = (ui, cropCoeff[i-1], refEvap[i-1], soilPars),
                              t_span=[0,samplingTimeInternal],y0=tuple(temp),method='BDF')
    temp = sol.y[:,-1]
    # ui = 0
headArray_1[i,:]=sol.y[:,-1]
volMoisture_1[i,:] = volMoistureAllNodes_1D(headArray_1[i,:], soilPars).ravel()

for i in range(2, totTimeSteps + 1):  # totTimeSteps + 1
    print(i)
    temp = headArray_1[i-1]
    ui = u_train[i-1]
    for t in range(int(samplingTime/samplingTimeInternal)): 
        sol = integrate.solve_ivp(ODE, args = (ui, cropCoeff[i-1], refEvap[i-1], soilPars),
                              t_span=[0,samplingTimeInternal],y0=tuple(temp),method='BDF')
        temp = sol.y[:,-1]
    # print(sol.y[:,-1])
        # ui = 0
    headArray_1[i,:]=sol.y[:,-1]
    volMoisture_1[i,:] = volMoistureAllNodes_1D(headArray_1[i,:], soilPars).ravel()

plt.figure(2)
plt.plot(volMoisture_1[2000:8000,19])
# np.savetxt("u_20hrLF_SRNN_100000_012_038.txt", u_train)
# np.savetxt("x_20hrLF_SRNN_100000_012_038.txt", headArray_1)
# np.savetxt("y_20hrLF_SRNN_100000_012_038.txt", volMoisture_1)

# np.savetxt("u_vali_80hrLF_2layered_1000_010_034.txt", u_train)
# np.savetxt("x_vali_50hrLF_2layered_500_011_035.txt", headArray_1)
# np.savetxt("y_vali_80hrLF_2layered_1000_010_034.txt", volMoisture_1)
# %% 
past = 20
future = 1
y_index = 19

def multivariate_data(dataset, target, start_index, end_index, history_size,
                      target_size,u_shape): #, single_step=False):
    data = []
    labels = []

    start_index = start_index + history_size
    if end_index is None:
        end_index = len(dataset) - target_size
        
    for i in range(start_index, end_index):
        temp_p = dataset[i-history_size:i,:]
        for j in range(target_size-1):
        # temp_p = np.concatenate(((dataset[i-history_size:i,:]).flatten(), 
        #                            (dataset[i:i+target_size-1,:u_shape]).flatten()))
            temp = np.array([dataset[i+j,0], dataset[i,1]]).reshape((1,2))
            temp_p = np.concatenate((temp_p,temp))
        data.append(temp_p)
        labels.append((target[i:i+target_size].flatten()))
    

    return np.array(data), np.array(labels)
# y = volMoisture_1[:,19]
# u = u_train
u = np.loadtxt("Modeling/u_20hrLF_SRNN_100000_012_038.txt")[:]
y = np.loadtxt("Modeling/y_20hrLF_SRNN_100000_012_038.txt")[:,y_index]
uv = np.loadtxt("u_vali_80hrLF_2layered_1000_010_034.txt")
yv = np.loadtxt("y_vali_80hrLF_2layered_1000_010_034.txt")[:,y_index]

# uv = np.loadtxt('u_NN_train_crop_2hr_impulse_f30.txt')[134000:-1]
# yv = np.loadtxt("y_NN_train_crop_2hr_impulse_f30.txt")[134001: ,y_index]
np.random.seed(5)
yn_train = 0.1*(max(y) - min(y))*np.random.choice([-1,1],(y.shape))*np.random.random((y.shape))
y_n = y*(1+yn_train)

y_train = y[1:]
x_train = np.zeros((u.shape[0],2))
x_train[:,0] = u
x_train[:,1] = y_n[:-1]

scaler_x = MinMaxScaler()
x_train = scaler_x.fit_transform(x_train)
scaler_y = MinMaxScaler()
y_train = scaler_y.fit_transform(y_train.reshape((y_train.shape[0],1)))
x_data, y_data = multivariate_data(x_train, y_train, 0, None, past, future,1)


yn = 0.2*(max(yv) - min(yv))*np.random.choice([-1,1],(yv.shape))*np.random.random((yv.shape))
yv_n = yv*(1+yn)

y_test = yv[1:]
x_test = np.zeros((uv.shape[0],2))
x_test[:,0] = uv
x_test[:,1] = yv_n[:-1]


scaler_xt = MinMaxScaler()
x_test = scaler_xt.fit_transform(x_test)
y_test = y_test.reshape((y_test.shape[0],1))
scaler_yt = MinMaxScaler()
y_test = scaler_yt.fit_transform(y_test)

x_t, y_t = multivariate_data(x_test, y_test, 0, None, past, future,1)

OUTPUT_SIZE = future*Ny

# %% Model 1 - First 120000 data points
model_1 = tf.keras.Sequential() 
model_1.add(tf.keras.layers.LSTM(5, input_shape=(past+future-1,2), activation = "sigmoid", return_sequences=True))
# model_1.add(tf.keras.layers.LSTM(3, activation = "sigmoid", return_sequences=True))
# model_1.add(tf.keras.layers.LSTM(5, activation = "sigmoid"))
# model_1.add(tf.keras.layers.Dense(5, input_shape=(past*2+future-1,), activation = "sigmoid"))  #, activation = 'sigmoid'
# model.add(Dense(50))
# model.add(tf.keras.layers.SimpleRNN(10, activation = 'tanh', return_sequences=True))
# model_1.add(tf.keras.layers.Dense(20, activation = 'sigmoid')) 
# model_1.add(tf.keras.layers.Dense(20, activation = 'sigmoid')) 
# model_1.add(tf.keras.layers.LSTM(20, return_sequences=True)) 
model_1.add(tf.keras.layers.LSTM(5, activation = "sigmoid"))  #, return_sequences=True
# model_1.add(tf.keras.layers.Dense(20, activation = 'tanh'))  # , activation = 'sigmoid',activation="selu"
# model.add(tf.keras.layers.Dense(1, activation = 'selu'))
# add output layer
# model_2.add(tf.keras.layers.Dense(5, activation = 'tanh')) 
model_1.add(tf.keras.layers.Dense(OUTPUT_SIZE, activation = "tanh"))

opt = tf.keras.optimizers.Adam(learning_rate=0.02)
model_1.compile(optimizer=opt, loss='mse', metrics=['mean_squared_error'])

epoch = 500

callback = tf.keras.callbacks.EarlyStopping(monitor='loss', patience=3)
Dense_history = model_1.fit(x_data,y_data.reshape((y_data.shape[0],y_data.shape[1])), epochs=epoch,
                               batch_size=200,
                                validation_split=0.15,
                           # steps_per_epoch = steps_per_epoch,
                                # validation_data = (x_t, y_t),
                                # validation_steps = validation_steps,
                               callbacks = [callback])


# %%
y_pt_1 = model_1.predict(x_t)

ypt_1 = scaler_yt.inverse_transform(y_pt_1)
# ypt_2 = scaler_yt.inverse_transform(y_pt_2)
# ypt = (ypt_1 + ypt_2)/2
y_act = scaler_yt.inverse_transform(y_t.reshape((y_t.shape[0], y_t.shape[1])))


plt.figure(3)
plt.plot(y_act[:,0], label = 'Actual Data')
plt.plot(ypt_1[:,0], label = 'Single-step-ahead Prediction')
# plt.plot(ypt_1[-500:,-1], label = 'Multi-step-ahead NN')
# plt.plot(ypt_2[30:,0], label = 'Multi-step-ahead NN -1st')
# plt.plot(ypt_2[:,-1], label = 'Multi-step-ahead NN - 30th')
plt.xlabel('Time Steps', **csfont, fontsize = 17)
plt.ylabel('Soil Moisture ($m^3/m^3$)', **csfont, fontsize = 17)
plt.legend(frameon=False, ncol=1, prop=font)

error = (sum((ypt_1 - y_act)**2/y_act**2))/ypt_1.shape[0]
print(error)
# %%
def sigmoid(x):
    return 1/(1+np.exp(-x))

def LSTMlayer_sig(warr, uarr, barr, x_t,h_tm1,c_tm1):
    
    '''
    c_tm1 = np.array([0,0]).reshape(1,2)
    h_tm1 = np.array([0,0]).reshape(1,2)
    x_t   = np.array([1]).reshape(1,1)
    
    warr.shape = (nfeature,hunits*4)
    uarr.shape = (hunits,hunits*4)
    barr.shape = (hunits*4,)
    
    '''
    
    s_t = (mtimes(x_t,warr) + mtimes(h_tm1,uarr) + barr.reshape((1, barr.shape[0])))
    hunit = uarr.shape[0]
    i  = sigmoid(s_t[:,:hunit])
    f  = sigmoid(s_t[:,1*hunit:2*hunit])
    # _c = np.tanh(s_t[:,2*hunit:3*hunit])
    _c = sigmoid(s_t[:,2*hunit:3*hunit])
    o  = sigmoid(s_t[:,3*hunit:])
    c_t = i*_c + f*c_tm1
    # h_t = o*np.tanh(c_t)
    h_t = o*sigmoid(c_t)
    return h_t, c_t

def LSTMlayer_tanh(warr, uarr, barr, x_t,h_tm1,c_tm1):
    
    '''
    c_tm1 = np.array([0,0]).reshape(1,2)
    h_tm1 = np.array([0,0]).reshape(1,2)
    x_t   = np.array([1]).reshape(1,1)
    
    warr.shape = (nfeature,hunits*4)
    uarr.shape = (hunits,hunits*4)
    barr.shape = (hunits*4,)
    
    '''
    
    s_t = (mtimes(x_t,warr) + mtimes(h_tm1,uarr) + barr.reshape((1, barr.shape[0])))
    hunit = uarr.shape[0]
    i  = sigmoid(s_t[:,:hunit])
    f  = sigmoid(s_t[:,1*hunit:2*hunit])
    _c = np.tanh(s_t[:,2*hunit:3*hunit])
    # _c = sigmoid(s_t[:,2*hunit:3*hunit])
    o  = sigmoid(s_t[:,3*hunit:])
    c_t = i*_c + f*c_tm1
    h_t = o*np.tanh(c_t)
    # h_t = o*sigmoid(c_t)
    return h_t, c_t




# %% Define sub model 3
model_1 = tf.keras.models.load_model("Modeling/LSTM_5s_5s_t_10nt")
# model_3 = model_1
def M(x):
    weightLSTM_1 = model_1.layers[0].get_weights()
    warr1, uarr1, barr1 = weightLSTM_1

    weightLSTM_2 = model_1.layers[1].get_weights()
    warr2, uarr2, barr2 = weightLSTM_2
    
    weightD = model_1.layers[2].get_weights()
    kernel3, bias3 = weightD
    
    n1 = 5
    # Initialise the model with the hidden states
    h_t1 = SX.zeros((1, n1))
    # Long term memory
    C_t1 = SX.zeros((1, n1))
    for i in range(x.shape[0]):
        xi = x[i,:].reshape((1,2))
        L1 = LSTMlayer_sig(warr1, uarr1, barr1, xi, h_t1[i,:], C_t1[i,:])
        h_t1 = vertcat(h_t1, L1[0])
        C_t1 = vertcat(C_t1, L1[1])
    
    n2 = 5
    # Initialise the model with the hidden states
    h_t2 = SX.zeros((1, n2))
    # Long term memory
    C_t2 = SX.zeros((1, n2))
    for i in range(h_t1.shape[0]):
        xi = h_t1[i,:].reshape((1,n1))
        L2 = LSTMlayer_sig(warr2, uarr2, barr2, xi, h_t2[i,:], C_t2[i,:])
        h_t2 = vertcat(h_t2, L2[0])
        C_t2 = vertcat(C_t2, L2[1])
    #  array([-2.7764497, -2.6578658, -2.6829786, -2.7356846, -2.73696  ]
    # Dense layer
    out2 = []
    for i in range(bias3.shape[0]):
        temp2 = mtimes(h_t2[-1,:], kernel3)[i]
    temp2 = bias3.reshape((1, bias3.shape[0])) + mtimes(h_t2[-1,:], kernel3)
    out2 = tanh(temp2)     
    
    return out2
# %% Multi-step ahead prediction performance
a = 900
t_pred = 20  # 30 sampling times 
umax = max(uv)
umin = min(uv)
u_sc = (uv - umin)/(umax - umin)
error = np.zeros(a)
y_save = np.zeros((a,t_pred))

for t0 in range(43, a+43):
    print(t0)
    y_sim_1 = np.zeros((t_pred,1))
    NN0_1 = x_t[t0,:]
    # y_sub = np.array([DM(M1(NN0_1)), DM(M2(NN0_1)), DM(M1(NN0_1))])  #, DM(M4(NN0_1)), DM(M5(NN0_1))])
    y_sim_1[0,:] = np.array(DM(M(NN0_1))).flatten()
    u_sim = x_t[t0:t0+t_pred,-1,0]
    for i in range(t_pred-1):
        NN0_1 = np.concatenate((NN0_1[1:], np.array([u_sim[i], y_sim_1[i,0]]).reshape((1,2))))
        # y_sub = np.array([DM(M1(NN0_1)), DM(M2(NN0_1)), DM(M1(NN0_1))])  #, DM(M4(NN0_1)), DM(M5(NN0_1))])
        y_sim_1[i+1,:] = np.array(DM(M(NN0_1))).flatten()
    yp_1 = scaler_yt.inverse_transform(y_sim_1).flatten()
    y_save[t0-43,:] = yp_1
    error[t0-43] = sum((y_act[t0:t0+t_pred].flatten() - yp_1)**2/y_act[t0:t0+t_pred].flatten())
# yp_2 = scaler_yt.inverse_transform(y_sim_2)
# yp = (yp_1 + yp_2)/2

# %% Plot
# y_save = np.loadtxt("Modeling/y_pred_80hrLF_2layered_1000_010_034_20nv_LSTM.txt")
y_save = np.loadtxt("Modeling/y_pred_80hrLF_2layered_1000_010_034_20nv_2layer_sub_noise.txt")
font = font_manager.FontProperties(family='Times New Roman',
                                   # weight='bold',
                                    style='normal', size=18)
csfont = {'fontname':'Times New Roman'}

a=900
b = 330-43
c = 385-43
tplot = np.arange(a+19)
fig,ax1 = plt.subplots()
# ax1.figure(5)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
ax1.plot(tplot[:a], y_act[43:a+43], label = 'Actual Trajectory', color = "k", linewidth = 3)
ax1.plot(tplot[:a], y_save[:a,0], label = "Single-step-ahead") 
# ax1.plot(tplot[4:], y_save[:-4,4], label = "Five-step-ahead", color = "r")  
ax1.plot(tplot[9:a+9], y_save[:a,9], label = "Ten-step-ahead", color = "r", linestyle = "dashdot") 
ax1.plot(tplot[19:], y_save[:a,19], label = "Twenty-step-ahead", color = "y", linestyle = "dashed") 
ax1.legend(frameon=False, ncol=1, prop=font)
plt.xlabel("Time Steps", **csfont, fontsize = 17)
plt.ylabel('Soil Moisture ($m^3/m^3$)', **csfont, fontsize = 17)


ax2 = plt.axes([0,0,1,1])
ip = InsetPosition(ax1, [0.55,0.75,0.15,0.2])
ax2.set_axes_locator(ip)
mark_inset(ax1, ax2, loc1=2, loc2=4, fc="none", ec='0.6')
ax2.plot(tplot[b:c+20], y_act[b+43:c+20+43], label = 'Actual Trajectory', color = "k", linewidth = 3)
ax2.plot(tplot[b:c], y_save[b:c,0], label = "Single-step-ahead") 
ax2.plot(tplot[b+9:c+9], y_save[b:c,9], label = "Ten-step-ahead", color = "r", linestyle = "dashdot") 
ax2.plot(tplot[b+19:c+19], y_save[b:c,19], label = "Twenty-step-ahead", color = "y", linestyle = "dashed") 
# ax2.legend(loc=0)
# ax2.set_xticklabels(ax2.get_xticks(), backgroundcolor='w')
ax2.tick_params(axis='x', which='major', pad=8)
# plt.suptitle('Composite Plot with original subsystems(Ca)', fontsize=15)
plt.show()

# plt.figure(4)
fig,ax1 = plt.subplots()
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
ax1.plot(tplot, y_act[:a], label = 'Actual Trajectory', color = "k", linewidth = 3)
ax1.plot(tplot, y_save_2[:a,0], label = "Single-step-ahead") 
# ax1.plot(tplot[4:], y_save[:-4,4], label = "Five-step-ahead", color = "r")  
ax1.plot(tplot[9:], y_save_2[:a-9,9], label = "Ten-step-ahead", color = "r", linestyle = "dashdot") 
ax1.plot(tplot[19:], y_save_2[:a-19,19], label = "Twenty-step-ahead", color = "y", linestyle = "dashed") 
ax1.legend(frameon=False, ncol=1, prop=font)
plt.xlabel("Time Steps", **csfont, fontsize = 18)
plt.ylabel('Soil Moisture ($m^3/m^3$)', **csfont, fontsize = 18)

ax2 = plt.axes([0,0,1,1])
ip = InsetPosition(ax1, [0.6,0.75,0.15,0.2])
ax2.set_axes_locator(ip)
mark_inset(ax1, ax2, loc1=2, loc2=4, fc="none", ec='0.6')
ax2.plot(tplot[b:c+20], y_act[b:c+20], label = 'Actual Trajectory', color = "k", linewidth = 3)
ax2.plot(tplot[b:c], y_save_2[b:c,0], label = "Single-step-ahead") 
ax2.plot(tplot[b+9:c+9], y_save_2[b:c,9], label = "Ten-step-ahead", color = "r", linestyle = "dashdot") 
ax2.plot(tplot[b+19:c+19], y_save_2[b:c,19], label = "Twenty-step-ahead", color = "y", linestyle = "dashed") 
# ax2.legend(loc=0)
# ax2.set_xticklabels(ax2.get_xticks(), backgroundcolor='w')
ax2.tick_params(axis='x', which='major', pad=8)
# plt.suptitle('Composite Plot with original subsystems(Ca)', fontsize=15)
plt.show()
# np.savetxt("y_pred_80hrLF_2layered_1000_010_034_20nv_LSTM.txt",y_save)
# model_1.save("LSTM_5s_5s_t_10nt_2")


e_01 = np.sqrt(sum((y_act[:900,0] - y_save_1[:,0])**2)/900)/(max(y_act[:900,0]) - min(y_act[:900,0])) #(max((y_act[:900,:] - y_save[:,:]).flatten()))
e_51 = np.sqrt(sum((y_act[4:904,0] - y_save_1[:,4])**2)/900)/(max(y_act[:900,0]) - min(y_act[:900,0])) 
e_101 = np.sqrt(sum((y_act[9:909,0] - y_save_1[:,9])**2)/900)/(max(y_act[:900,0]) - min(y_act[:900,0]))
e_201 = np.sqrt(sum((y_act[19:919,0] - y_save_1[:,19])**2)/900)/(max(y_act[:900,0]) - min(y_act[:900,0]))

print(e_01*100)
print(e_51*100)
print(e_101*100)
print(e_201*100)

e_02 = np.sqrt(sum((y_act[:900,0] - y_save_2[:,0])**2)/900)/(max(y_act[:900,0]) - min(y_act[:900,0])) #(max((y_act[:900,:] - y_save[:,:]).flatten()))
e_52 = np.sqrt(sum((y_act[4:904,0] - y_save_2[:,4])**2)/900)/(max(y_act[:900,0]) - min(y_act[:900,0])) 
e_102 = np.sqrt(sum((y_act[9:909,0] - y_save_2[:,9])**2)/900)/(max(y_act[:900,0]) - min(y_act[:900,0]))
e_202 = np.sqrt(sum((y_act[19:919,0] - y_save_2[:,19])**2)/900)/(max(y_act[:900,0]) - min(y_act[:900,0]))

print(e_02*100)
print(e_52*100)
print(e_102*100)
print(e_202*100)

# e_0 = np.sqrt(sum((y_act[:900,0] - y_save_2[:,0])**2)/900)/(max(y_act[:900,0] - y_save_2[:,0])) # 0.0009256418306921692
# e_5 = np.sqrt(sum((y_act[4:904,0] - y_save_2[:,4])**2)/900)/(max(y_act[:900,0] - y_save_2[:,0]))  # 0.0016076390558023086
# e_10 = np.sqrt(sum((y_act[9:909,0] - y_save_2[:,9])**2)/900)/(max(y_act[:900,0] - y_save_2[:,0]))  # 0.0020653705982084693
# e_20 = np.sqrt(sum((y_act[19:919,0] - y_save_2[:,19])**2)/900)/(max(y_act[:900,0] - y_save_2[:,0]))  # 0.0022271377200664584

# %%
u = np.loadtxt("Modeling/u_20hrLF_SRNN_100000_012_038.txt")[:]
y = np.loadtxt("Modeling/y_20hrLF_SRNN_100000_012_038.txt")[:,y_index]
yn = 0.1*(max(y) - min(y))*np.random.choice([-1,1],(y.shape))*np.random.random((y.shape))
y_n = y*(1+yn)
fig, axs = plt.subplots(2,sharex = True, figsize=(11,9))

a = 200
b = 800
t = np.arange(b-a+1)
axs[0].step(t[:b-a],-u[a:b], color = "k") 
# axs[1].plot(t[:a],y[:a]) #, label = 'Product Concentration')
axs[0].set_ylabel('Irrigation Rate ($m/s$)', fontsize='xx-large', labelpad=20, fontproperties=font)
axs[1].set_xlabel("Time Steps", fontsize='xx-large', fontproperties=font)
axs[1].set_ylabel('Soil Moisture ($m^3/m^3$)', fontsize='xx-large', labelpad=5, fontproperties=font)
# for ax in axs.flatten():
#     labels = ax.get_xticklabels() + ax.get_yticklabels()
#     [label.set_fontname('Times New Roman') for label in labels]
#     [label.set_fontsize('xx-large') for label in labels]
    

csfont = {'fontname':'Times New Roman'}
# plt.figure(2,figsize=(6, 4.5))
# plt.style.use('grayscale')
axs[1].plot(t, y[a:b+1], label = "Nominal", color = "k",linewidth = 3)
axs[1].plot(t, y_n[a:b+1], label = "10% White Noise", color = "r")
# plt.xlabel('Time (min)', **csfont, fontsize = 17)
# plt.ylabel('Product concentration ($mole/m^3$)', **csfont, fontsize = 17)
font = font_manager.FontProperties(family='Times New Roman',
                                   # weight='bold',
                                    style='normal', size=17)
axs[1].legend(frameon=False, ncol=1, loc="upper right", prop=font)
plt.show()

# tp = 200
# plt.figure(5)
# plt.plot(np.arange(100), y_act[tp-50:tp+50])
# plt.plot(np.arange(50,70), y_save[tp,:])
# plt.plot(np.arange(51,71), y_save[tp+1,:])
# plt.plot(np.arange(52,72), y_save[tp+2,:])
# plt.plot(np.arange(53,73), y_save[tp+3,:])
